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Abstract:  We  propose  a  new  approach  to  robot  path  planning  which  consists  of  incrementally 
building  a  graph  connecting  the  local  minima  of  a  potential  function  defined  in  the  robot’s  configu¬ 
ration  space  and  concurrently  searching  this  graph  until  a  goal  configuration  is  attained.  Unlike  the 
so-called  “global”  path  planning  methods,  this  approach  does  not  require  an  expensive  computation 
step  before  the  search  for  a  path  can  actually  start.  On  the  other  hand,  it  searches  a  graph  that 
is  usually  much  smaller  than  the  graph  searched  by  the  so-called  “local”  methods.  We  describe  a 
collection  of  specific  techniques  that  allow  to  engineer  various  implementations  of  our  path  planning 
approach.  The  purpose  of  these  techniques  is  to  (1)  construct  “good”  potential  fields  and  (2)  ef¬ 
ficiently  escape  their  local  minima  (i.e.,  efficiently  build  the  local  minima  graph).  They  are  based 
on  the  use  of  multiscale  pyramids  of  bitmap  arrays  for  representing  both  the  robot’s  workspace  and 
configuration  space.  This  distributed  representation  makes  it  possible  to  construct  potential  fields 
numerically,  rather  than  analytically ,  in  relation  to  other  efficient  numerical  techniques.  We  have 
implemented  these  techniques  in  a  path  planner,  which  turns  out  to  be  both  very  fast  and  capable 
of  handling  many  degrees  of  freedom.  The  planner  has  solved  a  variety  of  problems.  Some  of  them 
are  far  beyond  the  capabilities  of  previously  existing  planners. 
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contract  N00014-88-K-0620  (Office  of  Naval  Research),  CIFE  (Center  for  Integrated  Facility  Engineering), 
CIS  (Center  of  Integrated  Systems),  and  DEC  (Digital  Equipment  Corporation). 


1  Introduction 

In  this  paper  we  describe  several  numerical  potential  field  techniques  for  robot  motion  plan¬ 
ning.  We  have  implemented  these  techniques  in  a  path  planner  which  turns  out  to  be  both 
very  fast  and  capable  of  handling  many  degrees  of  freedom  (DOF’s).  In  particular,  the 
planner  has  been  able  to  plan  the  motions  of  mobile  robots  with  3  DOF’s  (two  translations 
and  one  rotation)  two  orders  of  magnitude  faster  than  most  other  existing  planners.  It 
has  also  generated  paths  of  robots  with  many  DOF’s  in  reasonable  amount  of  time.  For 
example,  difficult  paths  for  a  non-serial  manipulator  with  10  joints  (some  revolute,  some 
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prismatic)  were  produced  in  1  to  5  minutes1.  Paths  were  also  generated  in  a  3D  workspace 
for  a  manipulator  with  31  DOF’s  in  about  15  minutes  of  computation  time.  These  results 
are  far  beyond  the  capabilities  of  previously  implemented  planners. 

The  problem  of  generating  collision-free  paths  has  attracted  considerable  interest  during  the 
past  years  [Schwartz,  Hopcroft  and  Sharir,  1987]  [Schwartz  and  Yap,  1987].  By  simplifying, 
we  can  say  that  two  extreme  approaches  have  been  explored,  the  “global”  one  and  the 
“local”  one.  The  global  approach  consists  of  first  constructing  a  concise  representation 
of  the  connectivity  of  the  set  of  collision-free  configurations  of  the  robot  in  the  form  of  a 
“connectivity  graph”  and  then  searching  this  graph  for  a  path.  Various  techniques  have 
been  devised,  e.g.:  exact  cell  decomposition  [Schwartz  and  Sharir,  1983],  approximate  cell 
decomposition  [Brooks  and  Lozano-Perez,  1983]  [Faverjon,  1986],  retraction  on  a  curve 
network  [O’Dunlaing,  Sharir  and  Yap,  1983].  The  local  approach  consists  of  searching  a 
grid  placed  onto  the  robot’s  configuration  space  [Donald,  1984].  Heuristics  computed  from 
partial  information  about  the  geometry  of  the  configuration  space  are  used  to  guide  the 
search  (and  eventually  to  construct  the  grid  itself).  A  widely  used  heuristic  technique 
guides  the  search  for  a  path  along  the  flow  of  the  negated  gradient  vector  field  generated 
by  an  artificial  potential  field  [Khatib,  1986]  [Faverjon  and  Tournassoud,  1987]. 

The  drawback  of  the  global  approach  is  that  it  requires  an  expensive  precomputation  step  - 
the  construction  of  the  connectivity  graph  —  before  the  search  for  a  path  can  actually  start. 
Since  the  time  required  by  the  computation  of  the  connectivity  graph  is  typically  exponential 
in  the  dimension  n  of  the  configuration  space  (i.e.,  the  number  of  DOF’s),  the  approach  is 
untractable  even  for  reasonably  small  values  of  n.  To  our  knowledge,  no  effective  planner 
has  been  implemented  using  this  approach  with  n  >  4.  Instead,  the  local  approach  requires 
no  expensive  precomputation  step  before  starting  the  search  of  a  path.  Consequently, 
in  favorable  cases,  it  runs  substantially  faster  than  any  global  method.  But,  since  the 
search  graph  (i.e.,  the  grid)  is  considerably  larger  than  the  connectivity  graph  searched  by 
global  methods,  it  may  require  much  more  time  than  global  methods  in  unfavorable  cases. 
In  order  to  deal  with  this  difficulty,  local  methods  need  powerful  heuristics  to  guide  the 
search.  But  known  such  heuristics  have  the  drawback  of  eventually  leading  the  search  to 
dead-ends,  for  instance  local  minima  of  the  potential  field,  which  may  be  quite  difficult 
to  escape.  One  may  think  of  constructing  a  potential  field  with  no  other  local  minimum 
than  the  goal  configuration  (in  the  connected  subset  of  the  free  space  containing  the  goal 
configuration),  but  the  analytical  definition  of  such  a  potential  turns  out  to  be  very  difficult 
(e.g.,  see  [Koditschek  and  Rimon,  1989]).  Futhermore,  even  if  a  definition  was  available, 
it  is  likely  that  its  computation  would  constitute  an  expensive  precomputation  step  before 
path  generation,  similar  in  drawback  to  the  construction  of  a  connectivity  graph. 

There  are  many  reasons  motivating  the  development  of  a  fast  path  planner  capable  of  dealing 
with  many  DOF’s.  For  instance,  a  (semi-)autonomous  robot  will  have  to  generate  its  paths 
on-line,  based  on  its  current  model  of  the  world,  and  to  react  rapidly  to  contingencies. 
Although  one  may  envision  a  robot  that  learns  about  its  workspace  and  memorizes  a  variety 

1Most  of  the  experiments  reported  in  this  paper  were  carried  out  on  a  DEC  3100  MIPS-based  workstation 
using  simulated  robots. 
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of  typical  paths,  it  is  by  far  more  appropriate  to  have  a  fast  “real-time”  path  planner2. 
Robots  will  typically  combine  one  or  several  arms  mounted  on  a  mobile  vehicle;  for  instance, 
space  robots  may  consist  of  several  manipulator  arms  attached  to  a  free-flying  platform 
[Schneider,  1989].  Planning  the  paths  of  such  robots  will  require  to  be  able  to  handle  many 
DOF's,  specially  if  cooperation  among  them  is  needed.  One  may  argue  that,  most  of  the 
time,  at  every  instant,  the  planner  has  to  worry  only  about  a  subset  of  these  DOF's.  But 
determining  which  DOF’s  are  important  at  every  instant  is  also  part  of  planning  and  hence 
should  be  handled  by  the  planner.  Fast  path  planning  capabilities  may  also  be  extremely 
useful  for  off-line  programming  of  industrial  and  construction  robots,  which  requires  quick 
man-machine  graphical  interaction. 

We  propose  a  new  approach  to  path  planning  that  attempts  to  combine  the  advantages  of 
both  the  local  approach  (avoid  expensive  precomputation)  and  the  global  approach  (search 
a  concise  graph).  This  approach  consists  of  incrementally  building  a  graph  connecting  the 
local  minima  of  a  potential  function  defined  in  the  configuration  space  of  the  robot  and 
concurrently  searching  this  graph  until  a  goal  configuration  is  attained.  The  graph  of  local 
minima  plays  a  role  similar  to  that  of  the  connectivity  graph  in  the  global  approach.  The 
major  difference,  however,  is  that  it  is  constructed  in  an  incremental  fashion  during  the 
search.  Hence,  our  approach  does  not  require  an  expensive  precomputation  step,  although 
it  definitely  searches  a  much  smaller  graph  than  the  discretization  grid  thrown  over  config¬ 
uration  space.  In  this  paper,  we  describe  a  collection  of  specific  techniques  that  allow  to 
engineer  various  implementations  of  this  path  planning  approach3.  The  purpose  of  these 
techniques  is  to  (1)  construct  “good”  potential  fields  and  (2)  efficiently  escape  their  local 
minima  (i.e.,  efficiently  build  the  local  minima  graph).  They  are  based  on  the  use  of  multi¬ 
scale  pyramids  of  bitmap  arrays  for  representing  both  the  workspace  and  the  configuration 
space  of  the  robot.  These  representations  allow  us  to  construct  potential  fields  numerically , 
rather  than  analytically ,  in  relation  to  other  efficient  numerical  techniques  (e.g.,  collision 
checking,  valley  tracking). 

In  Section  2,  we  describe  the  hierarchical  bitmap  representations  of  the  workspace  and 
the  configuration  space  of  a  robot.  In  Section  3,  we  propose  several  ways  of  constructing 
numerical  potential  fields  in  the  robot’s  configuration  space.  In  Sections  4  through  7,  we 
present  four  path  planning  techniques  -  respectively  called  “best-first  motion”,  “random 
motion”,  “valley-guided  motion”  and  “constrained  motion”  techniques  -  which  are  based  on 
different  ways  of  escaping  the  local  minima  of  these  potential  fields.  We  have  implemented 
all  these  techniques  and,  for  each  one,  we  give  experimental  results.  Each  technique  admits 
many  straightforward  variants  and,  in  the  course  of  our  experimentation,  we  ran  several 
variants  successfully.  Hence,  within  some  limits,  the  algorithms  presented  in  this  paper 
may  be  adapted  so  that  they  better  fit  the  characteristics  of  a  specific  application  domain 
and  computing  system. 

2To  that  respect,  it  is  interesting  to  remember  that  not  so  many  years  ago,  it  was  proposed  to  compute 
the  inverse  kinematics  of  a  manipulator  by  storing  the  numerical  values  of  the  inverse  Jacobian  matrix  of  the 
forward  kinematic  map  at  many  configurations  of  the  manipulator.  Today,  the  computation  of  the  inverse 
kinematics  is  routinely  done  in  real-time  without  having  to  store  inverse  Jacobian  matrices. 

3 Some  of  these  techniques  were  previously  presented  in  [Barraquand,  Langlois  and  Latombe,  1989]  and 
[Barraquand  and  Latombe,  1989]. 
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2  Bitmap  Representations 

2.1  Basic  Terminology  and  Notations 

We  denote  by  A  the  robot  and  W  its  workspace.  A  standard  Cartesian  coordinate  frame, 
denoted  by  Tyy,  is  embedded  in  W. 

A  configuration  of  A  is  a  specification  of  the  position  of  every  points  in  A  with  respect 
to  jT-yy.  The  configuration  space  of  A  is  the  space,  denoted  by  C,  of  all  the  possible 
configurations  of  A.  The  subset  of  W  occupied  by  A  at  configuration  q  is  denoted  by  .4(q). 

The  workspace  contains  a  finite  number  of  obstacles  denoted  by  B{}  with  i  =  1, r.  We 
denote  by  YJ empty  the  region  W  -  Ui=i  ft-  Each  obstacle  Bi  maps  in  C  to  the  set  CBi  of 
configurations,  called  C-obstacle,  where  the  robot  and  the  obstacles  intersect,  i.e.: 

Cft  =  {qec/-4(q)nft  #0}. 

The  subset  of  configurations  where  the  robot  and  the  obstacles  have  no  intersection,  i.e.: 

Cjree  =  C  -  \J  CBt 

»=  1 

is  called  the  free  space.  A  collision-free  path  (more  simply,  a  free  path)  between  two 
configurations  qinit  and  q goal  is  any  continuous  map  r  :  [0, 1]  — »  C/ree>  such  that  r(0)  =  q init 
and  r(l)  =  q goal- 

The  configuration  space  C  is  a  nD  manifold  [Arnold,  1978].  For  example,  most  mobile  robots 
can  be  modelled  as  a  2D  object  A  that  can  both  translate  and  rotate  in  the  plane.  Then, 
C  is  the  3D  manifold  R2  x  51,  where  S1  is  the  unit  circle.  In  the  case  of  a  manipulator 
arm  with  n  revolute  joints,  C  is  the  nD  manifold  (S1)71,  or  a  subset  of  that  space  when  the 
motion  of  each  joint  is  limited  by  mechanical  stops. 

Throughout  this  paper,  we  represent  a  configuration  q  of  A  by  a  list  of  n  parameters, 
(<Zi ,...,gn),  where  n  is  the  dimension  of  C .  This  parameterization  is  eventually  augmented 
with  modular  arithmetic  for  some  of  the  angular  parameters.  Hence,  we  represent  C  as 
a  nD  Cartesian  space.  For  example,  C  =  R2  x  S1  may  be  represented  as  R2  x  H/2'jtZ. 
Any  configuration  q  is  then  parameterized  by  (x,y,  0),  where  (x,y)  €  R2  and  d  €  [0 , 2 7r) 
with  modulo  2tt  arithmetic.  Similarly,  C  =  (5x)n  may  be  represented  as  (R/27T Z)n .  Any 
configuration  is  parameterized  by  (gi,...,gn)  with  qi  £  [0, 2tt),  for  every  i  €  [l,n],  and 
modulo  2tv  arithmetic  on  every  g*. 

For  each  point  p  £  *4,  one  can  consider  the  geometrical  application  that  maps  any  configu¬ 
ration  q  =  (gi,  ...,gn)  G  C  to  the  position  x  G  W  of  p  in  the  workspace.  This  map: 

X  :  A  x  C  W 

Cp,  q)  ^  x  =  X(p,q) 


is  called  forward  kinematic  map. 
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Figure  1:  2562  Bitmap  Description  of  a  Workspace 


2.2  Workspace  Bitmap 

In  the  following,  we  make  the  reasonable  assumption  that  >V  is  a  bounded  subset  of  Rd, 
with  d  =  2  or  3.  W  is  modelled  as  a  multiscale  pyramid  of  dD  bitmap  arrays.  At  each 
resolution  level,  the  array  is  represented  by  a  function: 

BM  :  W  —  {1,0} 

x  (-»•  BM(x ) 

in  such  a  way  that  the  subset  of  points  x  such  that  BM[x)  =  1  represents  the  workspace 
obstacles  and  the  subset  of  points  x  such  that  BM(x)  =  0  represents  the  empty  part  of  the 
workspace,  i.e.  W empty 

The  distance  between  the  centers  of  two  consecutive  cells  in  the  same  line  or  the  same 
column  of  an  array  is  constant.  Hence,  at  any  level  of  resolution,  the  midpoints  of  the  cells 
of  the  bitmap  array  form  a  regular  grid  denoted  by  QW .  The  subset  of  the  grid  where 
BM  evaluates  to  0  is  denoted  by  GWempty  For  any  integer  k  €  [l,r],  the  £-heighborhood 
of  a  point  x  in  a  grid  of  dimension  r  is  defined  as  the  set  of  points  in  the  grid  having 
at  most  k  coordinates  differing  from  those  of  x,  the  amount  of  the  difference  (if  any) 
along  any  axis  being  the  discretization  step  along  that  axis.  In  GW,  we  always  use  the 
1-neighborhood,  except  noticed  otherwise  (e.g.,  to  track  a  curve  or  a  surface).  In  a  2D 
workspace  grid,  this  means  that  each  point  x  =  (i,j)  6  GW  has  a  maximum  of  4  neighbors: 

{(*  ~  -  l),(i,j  +  1)}.  In  a  3D  workspace  grid,  each  point  has  up  to  6 

neighbors. 

Most  of  the  experiments  reported  in  this  paper  were  carried  out  in  a  2D  workspace  (d  =  2). 
In  those  experiments,  the  coarsest  level  of  resolution  in  the  pyramid  was  162  and  the  finest 
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1.  For  every  x  €  QW  empty,  set  di(x)  to  infinity  (i.e.,  a  large  number  M). 

2.  Scan  QW  and  identify  every  point  x  such  that  BM(x)  =  1  and  one  of  its  neighbors  is 
in  GW  empty  Set  L0  to  the  list  of  these  points.  Include  the  points  forming  the  frame 
boundary  of  GW  in  Lq.  For  every  point  x  in  Lq,  set  di(x)  to  0.  Set  i  to  0. 

3.  For  i  =  0,1,2,...,  while  L,-  until  Li  is  empty,  do:  initialize  to  the  empty  list;  for 
every  point  x  in  £,♦,  for  every  neighbor  y  of  x  in  GW  empty,  if  di(y)  =  M  then  set  di(y) 
to  i  +  1  and  insert  y  at  the  end  of  L,-+i . 


Figure  2:  Algorithm  Computing  the  d\  Map 

5122.  The  bitmap  representation  of  a  particular  workspace  at  the  2562  resolution  is  shown 
in  Figure  1  (1  =  black;  0  =  white).  The  workspace  representation  is  given  to  the  planner 
at  the  finest  level  of  resolution.  The  other  levels  are  automatically  derived  from  it  in  a 
conservative  fashion,  so  that  {x  £  W  /  BM(x)  =  0}  C  Wempty  The  scaling  factor  between 
two  successive  levels  of  resolution  is  2,  but  a  different  factor  could  have  been  chosen. 

In  preparation  to  other  algorithms  to  be  described  later,  the  planner  computes  the  dis¬ 
cretized  L 1  distance  c?i(x)  from  every  point  x  €  GW  empty  to  the  obstacles.  This  computa¬ 
tion  is  performed  according  to  the  following  “wavefront  expansion”  algorithm4.  First,  the 
points  in  the  boundary  of  the  obstacles  are  identified  and  the  value  of  d\  at  these  points  is 
set  to  zero  (the  points  in  the  contour  of  the  bitmap  are  also  included  as  boundary  points). 
Next,  the  value  of  d\  is  set  to  1  at  all  the  neighbors  of  the  boundary  points  in  QW empty]  to  2 
at  the  neighbors  of  these  new  points  (if  not  yet  computed);  etc  ...  The  algorithm  terminates 
when  all  QW  empty  has  been  explored.  A  more  formal  description  of  the  algorithm  is  given 
in  Figure  2. 

The  time  complexity  of  this  algorithm  is  linear  in  the  number  of  points  in  the  grid  QW. 
Figure  3  displays  contours  of  d\  for  the  workspace  of  Figure  1. 

2*3  Configuration  Space  Bitmap 

Since  we  discretize  the  workspace  in  a  hierarchical  fashion,  it  is  consistent  to  discretize  the 
configuration  space  C  as  well,  by  constructing  another  multiresolution  grid  pyramid.  This 
pyramid  has  as  many  levels  of  resolution  as  the  workspace  pyramid  and  the  resolutions  at 
each  level  of  the  two  pyramids  are  tightly  related,  as  developed  below.  At  each  level  of 
resolution,  we  denote  by  QC  the  grid  representing  the  configuration  space  and  by  QCfree 
the  subset  of  the  grid  lying  in  the  free  space  C/ree. 

Let  us  denote  by  8  the  distance  between  two  adjacent  points  in  a  workspace  grid  QW .  In 
the  workspace  pyramid,  8  varies  between  8min  and  8rnax>  For  example,  let  us  assume  that 
we  have  a  workspace  represented  by  a  pyramid  of  arrays  whose  sizes  are  ranging  between 

4In  this  algorithm  we  normalize  the  distance  between  two  neighbors  in  the  grid  to  1. 
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Figure  3:  Contours  of  d\  in  a  Workspace 

162  and  5122.  If  the  distance  is  measured  in  percentage  of  the  workspace  diameter,  we  have 
^min  —  1/512  and  6 max  —  1/16. 

By  definition,  the  resolution  of  a  grid  is  the  logarithm  of  the  inverse  of  the  distance  between 
two  discretization  points,  in  the  base  defined  by  the  scaling  factor  between  two  successive 
resolution  levels  (2  in  our  implementation).  In  our  example,  the  resolution  r  hence  varies 
between  rmtn  =  —  log2(<5max)  =  4  and  rmax  ~  log2(<5mtn)  =  9. 

Remember  that  we  represent  the  configuration  space  C  as  a  subset  of  a  nD  Cartesian 
space  with  q  =  (gi,...,gn).  For  any  given  workspace  resolution,  say  r  =  -log2(<5),  the 
corresponding  resolution  Ri  =  -log2(Ai)  of  the  discretization  along  the  qi  axis  of  C  is 
chosen  in  such  a  way  that  a  modification  of  qi  by  A;  =  2~Ri  generates  a  “small  motion”  of 
the  robot  in  the  workspace.  By  “small  motion”,  we  mean  that  any  point  p  £  A  moves  by 
less  than  nbtol  x  <5,  where  nbtol  is  a  small  number  (typically,  1  or  2). 

The  relation  between  positions  of  robot  points  in  the  workspace  and  robot  configurations 
is  given  by  the  forward  kinematic  map  X(p ,  q).  For  every  point  p  £  A,  a  modification  of  qi 
by  A i  results  in  a  modification  of  each  coordinate  Xj  of  p  (j  £  [1,  d])  by: 

|f(P,q)A„ 

If  we  impose  all  workspace  motions  to  be  less  than  nbtol  x  8,  we  must  have: 

A,-  =  nbtol  x  8/  sup  ( ^-(p,  q)]  =  nbtol  x  8/J j 

P6^,qgc,je[i,d]  V  clqi  / 
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For  a  given  robot,  the  numbers  JgUp  are  generally  straightforward  to  compute.  This  leads 
to  compute  the  resolution  Ri  as: 

Ri  =  r  +  log2(JsUp)  -  \og2(nbtol) 

As  an  example,  let  us  consider  a  bar  of  length  L  moving  freely  in  a  2D  workspace.  We  can 
represent  a  configuration  of  the  bar  by  (xg,2/g,0),  where  xq  and  yo  are  the  coordinates  of 
the  center  of  gravity  and  6  is  the  orientation  of  the  bar.  Let  us  normalize  xg,  Vg ,  and  6  so 
that  their  values  range  between  0  and  1.  We  have 


TXG  —  JVG  —  1 
''sup  ‘'sup  A 


and: 

jLp  =  *L- 


If  we  set  nbtol  —  2,  we  get: 


RXq  —  RVg  —  r  —  1 


and: 


Re  =  r  +  log2(7ri)  -  1. 


This  means  that  we  need  21  =  2  times  less  samples  for  xq  and  yc  than  for  the  workspace 
representation  at  each  level  of  resolution,  and  2/{rcL)  less  samples  for  6. 


The  planning  techniques  described  below  do  not  construct  and/or  use  a  complete  represen¬ 
tation  of  the  configuration  space  grids,  nor  of  the  C-obstacles,  since  this  would  be  too  time 
and  space  consumming.  Indeed,  while  the  configuration  space  grids  implicitly  define  the 
search  space  of  the  planner,  in  most  practical  cases,  only  a  very  small  subset  of  the  grids 
will  be  explored.  It  is  clear,  however,  that  in  the  worst  case  the  grids  would  have  to  be 
exhaustively  explored,  requiring  exponential  time  computation  as  any  other  existing  path 
planning  method. 


3  Potential  Field  Construction 

3.1  Overview 

We  describe  below  several  ways  of  constructing  numerical  potential  fields  in  the  configu¬ 
ration  space  of  a  robot.  In  all  cases,  the  construction  consists  of  two  major  computing 
steps.  First,  potential  fields,  called  W-potentials,  are  computed  in  the  robot  workspace 
W.  Each  W-potential  applies  to  a  selected  point  in  the  robot,  called  control  point,  and 
pulls  this  point  toward  its  goal  position.  The  W-potentials  at  the  various  control  points  are 
then  combined  into  another  function,  called  C-potential,  which  is  defined  over  the  robot’s 
configuration  space. 

More  formally,  let  pi ,  i  =  1,  denote  the  control  points  in  the  robot.  Each  W-potential 
is  a  function: 

LTp^  .  X  £  W empty  LJpt-  0  IT 

The  C-potential  is  defined  as: 

U(q)  =  G(UP1  (X(Pl ,  q)), .  . . ,  Up,  (X(ps,  q))) 
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1.  For  every  x  G  GW  empty,  set  U  p(x)  to  infinity  (i.e.,  a  large  number  M). 

2.  Set  U p(xgoai)  to  0,  the  list  Lq  to  (x^oa/),  and  the  counter  i  to  0. 

3.  For  i  —  0, 1, ...,  until  Li  is  empty,  do:  initialize  L,-+i  to  the  empty  list;  for  every  point  x 
in  Li ,  for  every  neighbor  y  of  x  in  GW  empty,  if  Up(y)  =  M  then  set  Up(y)  to  i  +  1  and 
insert  y  at  the  end  of  Lt+i. 


Figure  4:  Algorithm  Computing  the  Simple  W-Potential 
where  G  is  called  the  arbitration  function. 

The  rationale  behind  this  two-step  computation  is  to  use  the  configuration  space  as  a  source 
of  inspiration  for  constructing  a  “good”  C-potential.  As  a  matter  of  fact,  the  W-potentials 
computed  by  the  algorithms  of  the  next  subsection  are  free  of  local  minima.  This  allows 
us  to  construct  C-potentials  which  avoid  the  robot  to  get  trapped  in  simple  concavities 
formed  by  the  obstacles.  However,  since  several  W-potentials  are  typically  combined  in 
order  to  construct  a  C-potential,  the  combination  may  still  have  spurious  local  minima. 
These  minima  result  from  the  fact  that  the  various  control  points,  which  are  related  by 
fixed  or  variable  kinematic  constraints,  are  concurrently  attracted  toward  their  respective 
goal  positions.  Hence,  they  are  competing  among  themselves  to  attain  their  respective 
goal  positions.  It  is  this  competition  which  may  create  unwanted  minima.  The  role  of  the 
function  G  is  to  arbitrate  in  this  competition,  hence  the  name  of  the  function.  Various 
arbitration  functions  are  possible,  resulting  in  more  or  less  numerous,  more  or  less  deep 
local  minima  (Subsection  3.3). 

3.2  Computation  of  W-Potentials 

We  propose  two  techniques  for  computing  local-minima-free  W-potentials.  Other  similar 
techniques  can  easily  be  developed.  In  the  rest  of  the  paper,  we  will  refer  to  the  two 
W-potentials  defined  below  as  the  “simple  W-potential”  and  the  “improved  W-potentiaP, 
respectively. 

3.2.1  Simple  W-Potential 

Let  p  be  a  control  point  in  the  robot  and  Up(x)  the  corresponding  W-potential.  We  want 
Up  to  be  free  of  local  minima,  i.e.  to  have  a  single  minimum  at  the  goal  position  of  p  in 
the  connected  subset  of  W/ree  containing  this  goal  position.  Indeed,  as  mentioned  above, 
we  think  that  this  is  a  major  heuristic  step  toward  the  construction  of  a  C-potential  with 
few  or  small  spurious  local  minima. 

Let  Xgoai  be  the  goal  position  of  p  in  W .  The  simple  W-potential5  Up  is  computed  as  follows. 
First,  the  value  of  Up  is  set  to  0  at  xgoai.  Next,  the  value  of  Up  is  set  to  1  at  the  neighbors 
of  x goal  in  Gempty]  to  2  at  the  neighbors  of  these  new  points  (if  not  yet  computed),  etc...  The 

5 A  similar  potential  has  previously  been  proposed  in  [Jarvis  and  Byrne,  1987]. 
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Figure  5:  Contours  of  the  Simple  W-Potential 


algorithm  terminates  when  all  GW  empty  has  been  explored.  A  more  formal  description  of 
the  algorithm  is  given  in  Figure  4.  At  every  point  x  E  GW  empty,  the  resulting  W-potential 
is  equal  to  the  L 1  length  of  the  minimal  length  path  connecting  x  to  x^oa/  through  GW  empty 
It  has  no  other  local  minimum  than  x.goai. 

This  algorithm  is  similar  in  structure  to  the  algorithm  computing  d\.  Its  time  complexity 
is  also  linear  in  the  number  of  points  of  GW  and  constant  in  the  number  and  shape  of  the 
obstacles.  Figure  5  displays  contours  of  Up  for  the  2D  workspace  of  Figure  1,  with  xgoal 
located  in  the  upper  lefthand  corner  of  the  bitmap  array.  The  computation  was  done  in  a 
fraction  of  a  second. 

A  property  of  the  Up  function  computed  as  above  is  the  following.  By  tracking  the  flow  of 
the  negated  gradient  vector  field  —  VUp  from  any  initial  point  Xim*,  we  obtain  a  path  that 
connects  xt*nt*t  to  x50a/.  In  addition,  this  path  is  the  shortest  path  for  the  L 1  distance  (at  the 
resolution  of  the  bitmap  array).  In  a  3D  workspace,  this  computation  may  be  preferable  to 
exact  methods,  since  the  problem  of  computing  the  exact  shortest  distance  in  a  polyhedral 
space  is  known  to  be  NP-hard  in  the  number  of  vertices  under  any  Lp  metric  [Canny,  1987]. 

Notice  that  the  algorithm  given  above  actually  computes  Up  only  in  the  connected  subset 
of  GW  empty  that  contains  the  goal  position  x50a/.  Hence,  if  the  initial  position  x2m*  of  p  is 
such  that  Up(x;m**)  =  M,  we  can  immediately  return  that  there  is  no  path  connecting  x;n^ 
to  'x.goal •  Step  2  of  the  above  algorithm  can  easily  be  modified  to  accommodate  the  case 
where  the  goal  of  p  is  a  subset  of  W. 


Figure  6:  Examples  of  Workspace  Skeletons 

3.2.2  Improved  W-Potential 

A  significant  drawback  of  the  simple  W-potential  is  that  it  induces  paths  which  typically 
graze  obstacles  in  the  workspace.  Since  several  W-potentials  usually  have  to  be  combined 
into  a  C-potential  function  U  attracting  the  robot  toward  a  goal  configuration,  the  indi¬ 
vidual  W-potentials  may  compete  in  such  a  way  that  they  produce  local  minima  of  U.  In 
order  to  reduce  the  risk  of  creating  local  minima  and  to  enlarge  the  maneuvering  space  of 
the  robot,  so  that,  if  local  minima  are  created,  they  can  be  more  easily  escaped,  we  propose 
an  improved  W-potential.  This  potential,  like  the  previous  one  is  free  of  local  minima.  In 
addition,  its  negated  gradient  pulls  the  control  point  p  toward  its  goal  position  along  a 
path  of  non-minimal  length,  which  stays  as  far  away  as  possible  from  the  obstacles.  This 
improved  W-potential  is  computed  in  three  stages. 

First,  the  discrete  Ll  distance  di(x)  from  every  point  x  £  QW  empty  to  the  obstacles  is 
computed  and  a  (d—  1)D  subset  7Z  of  QW  empty  from  the  map  d\  is  concurrently  extracted. 
This  subset  is  called  the  workspace  skeleton.  Second,  the  function  Up  in  the  skeleton  7 Z 
is  computed.  Third,  Up  is  computed  in  the  rest  of  QYd empty 

The  distance  di(x)  is  computed  using  the  algorithm  given  in  Figure  2.  The  workspace 
skeleton  71  is  extracted  during  the  computation  of  di  as  the  set  of  points  where  the  “waves” 
-  the  lists  Li  in  the  algorithm  -  issued  from  the  boundary  points  of  QV\> empty  meet.  This  is 
done  by  propagating  not  only  the  values  of  dx ,  but  also  the  points  in  the  discretized  boundary 
of  QYJ empty  that  are  at  the  origin  of  the  propagation.  For  each  point  x  in  QYdempty,  let 
O(x)  denote  this  point.  Notice  that  di(x)  is  the  Ll  distance  from  x  to  O(x).  At  step  2, 
we  set  O(x)  to  x  for  every  x  in  Lq.  At  step  3,  let  x  be  the  point  currently  considered  in  Lt- 
and  y  one  of  its  neighbors.  If  di(y)  =  M,  we  simply  set  O(y)  to  O(x);  otherwise,  if  the  L1 
distance  between  O(y)  and  O(x)  is  greater  than  some  threshold  (typically,  2),  we  include 
y  in  71. 

Figure  6  displays  the  skeletons  computed  for  several  2D  workspaces,  including  the  workspace 
shown  in  Figure  1.  Each  of  these  skeletons  is  a  kind  of  generalized  Voronoi  diagram  of 
y^empty  [Lee  and  Drysdale,  1981]  for  the  L1  distance.  It  is  also  similar  to  the  skeleton  ex- 
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lc  ( Initialization .) 

For  every  x  €  GW  empty,  set  Up(x)  to  infinity  (i.e.,  a  large  number  M).  Set  S  to  %. 

2.  ( Computation  of  the  augmented  skeleton  S  by  connecting  the  goal  to  the  skeleton .) 

Set  x  to  Xg0at .  While  x  £  11,  select  a  neighbor  y  of  x  having  the  largest  value  of  di, 
include  y  in  S  and  set  x  to  y, 

3.  ( Computation  of  the  W-potential  in  the  augmented  skeleton.) 

Set  Up(x^oa/)  to  0,  L  to  (xg0ai)  and  Sf  to  the  empty  list.  ( L  is  a  list  of  points  sorted 
by  decreasing  values  of  d\.)  Until  L  is  empty,  do:  remove  the  first  element  x  of  L  and 
insert  it  at  the  tail  of  S';  for  every  d-neighbor  y  of  x  in  S,  if  Up(y)  =  M  then  set  Up(y) 
to  Up(x)  +  1  and  insert  y  in  L.  (At  the  end  of  this  step,  S'  contains  all  the  points  in  S 
sorted  by  increasing  values  of  Up.) 

4.  ( Computation  of  the  W-potential  in  the  rest  ofGWernpty) 

Set  Lq  to  S'.  For  i  =  0, 1,2, ...,  until  L*  is  empty,  do:  initialize  L;+i  to  the  empty  list; 
for  every  point  x  in  for  every  neighbor  y  of  x  in  GWempty,  if  Up(y)  —  M  then  set 
Up(y)  to  Up(x)  +  1  and  insert  y  at  the  end  of  1. 


Figure  7:  Algorithm  Computing  the  Improved  W-Potential 

tracted  from  a  region  in  a  digitized  image  using  techniques  from  Mathematical  Morphology 
[Serra,  1982]. 

The  improved  W-potential  in  is  computed  as  follows  (see  Steps  2  and  3  in  Figure  7). 
First,  the  goal  position  xgoai  is  connected  to  71  by  a  path  M  following  the  gradient  of  dj. 
The  set  S  =  72.  U  M,  called  the  “augmented  skeleton”,  is  connected.  Then,  all  the  points 
in  S  are  labelled,  starting  at  xgoai ■  The  label  0  is  assigned  to  xgoai\  the  label  1  is  assigned 
to  every  d-neighbor6  of  xgoa\  in  S  which  is  at  maximum  distance  d\  of  the  obstacles;  the 
label  2  is  assigned  to  every  d-neighbor  of  a  previously  labelled  point  that  is  at  maximum 
distance  d\  of  the  obstacles;  etc  ...  The  algorithm  terminates  when  all  the  points  in  the 
augmented  skeleton  have  been  labelled.  At  each  iteration,  the  set  of  labelled  points  of 
5,  whose  neighbors  have  not  been  analyzed  yet,  is  sorted  by  decreasing  values  of  d\  in  a 
height-balanced  tree  [Aho,  Hopcroft  and  Ullman,  1983],  so  that  both  the  insertion  of  a  new 
point  and  the  removal  of  a  point  of  maximum  value  of  d\  are  done  in  logarithmic  time. 

The  W-potential  at  all  the  other  points  in  GW  empty  is  computed  using  a  wavefront  expansion 
algorithm  (see  Step  4  in  Figure  7)  starting  from  <S,  similar  to  the  algorithms  computing  the 
dx  map  and  the  simple  W-potential.  The  W-potential  at  each  1-neighbor  y  of  every  point  x 
in  S  is  set  to  Up(x)  +  1.  The  W-potential  at  the  neighbors  of  these  neighbors  is  iteratively 
incremented  until  all  the  points  in  QW empty  have  been  explored. 

Figure  8  shows  the  equipotential  contours  of  the  resulting  W-potential  for  the  workspace 
6 In  order  to  track  the  augmented  skeleton  reliably,  it  is  more  appropriate  to  use  the  d-neighborhood. 
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Figure  8:  Contours  of  Improved  W- Potential 


of  Figure  1.  This  W-potential  generates  no  stable  equilibrium  state.  Following  its  negated 
gradient  from  any  initial  position  xtnti  produces  a  path  that  first  connects  xtmt  to  'll,  then 
stays  as  far  as  possible  from  the  obstacles  by  following  the  safest  curve  of  H,  and  finally 
connects  'll  to  xgoai.  The  above  algorithm  only  computes  Up  over  the  connected  subset  of 
GW  empty  that  contains  xg0al-  It  can  be  adapted  to  the  case  where  the  goal  of  p  is  a  region 
in  W. 

The  complexity  of  computing  the  improved  W-potential  is  slightly  higher  than  for  the 
simple  one.  Let  a  be  the  number  of  points  in  the  bitmap  array,  and  b  the  number  of 
points  in  the  augmented  skeleton  <5.  The  complexity  of  the  algorithm  is  0(a  +  6  log  6).  For 
reasonable  workspaces,  however,  we  have  b  oc  a~  since  the  skeleton  is  a  (d-l)D  subset  of 
the  workspace.  For  such  workspaces,  the  complexity  is  0(a  +  a~  log  a),  hence  is  linear  in 
a.  For  the  2562  workspace  of  Figure  8,  the  computation  took  about  2  seconds,  including 
the  computation  of  d\  and  H. 

Variants  of  the  above  W-potential  are  easy  to  define  and  compute.  For  example,  in  step  4 
of  the  algorithm,  we  could  compute  Up(y)  =  Up(x)  +  l/di(x)  and  obtain  a  W-potential 
that  becomes  infinite  on  the  boundary  of  GWemVty  We  will  not  detail  the  cosmetics  of  the 
computation  of  numerical  W-potentials  further.  Our  point  is  simply  to  show  that  a  large 
family  of  W-potentials  with  various  properties  can  be  built  within  the  bitmap  numerical 
framework. 


13 


3.2.3  Remarks 


The  computation  of  di  (if  the  simple  W-potential  is  used)  and  d\  and  7 t  (if  the  improved 
W-potential  is  used)  is  not  local  and  therefore  must  be  done  prior  to  the  execution  of  the 
rest  of  any  planning  algorithm  based  on  our  path  planning  approach.  However,  the  time 
complexity  of  the  algorithm  is  only  a  function  of  the  number  of  points  in  QW.  It  is  constant 
in  both  the  number  and  the  shape  of  obstacles,  and  the  number  of  DOF’s  of  the  robot.  Its 
implementation  turns  out  to  be  quite  fast. 

From  a  conceptual  point  of  view,  neither  the  choice  of  the  L 1  metric,  nor  the  precise 
definition  of  the  skeleton  It  is  fundamental  to  the  rest  of  our  path  planning  approach. 
We  could  have  used  the  more  classical  L 2  distance  instead  and  computed  the  generalized 
Voronoi  diagram  for  that  metric.  However,  as  we  just  noticed,  the  construction  of  the 
workspace  potentials  over  QWempty  is  a  necessary  precomputation  before  the  rest  of  our 
path  planning  approach  can  be  executed.  The  computation  of  the  L 1  distance  is  faster  than 
the  computation  of  the  L 2  distance.  Notice  also  that,  unlike  the  above  computation  of  It, 
the  time  complexity  of  constructing  the  algebraic  description  of  the  L 2  Voronoi  diagram  of 
a  polygonal  workspace  increases  with  the  number  of  vertices  of  the  obstacles  [Leven  and 
Sharir,  1987].  We  experimented  with  variants  of  the  improved  W-potential  using  slightly 
different  constructions  of  the  workspace  skeleton  It.  In  general,  these  variants  produced 
similarly  satisfactory  results. 

3.3  Computation  of  C-Potentials 

The  C-potential  U  is  computed  as  a  combination: 

U(q)  =  G(UP1  (*0* ,  q)), . . . ,  UPs  (X(p, ,  q))) 

of  the  W-potentials  Up. ,  i  =  1, ...,  s,  defined  for  several  control  points  pi.  This  combination 
concurrently  attracts  the  different  points  pi  towards  their  respective  goal  positions.  U(q) 
attains  its  minimal  value,  i.e.  0,  when  all  the  control  points  are  at  their  goal  positions.  At 
any  other  configuration,  it  is  strictly  positive. 

The  control  points  are  those  used  to  input  the  description  of  the  goal  configuration  of  the 
robot.  (By  definition,  the  robot  is  at  a  goal  configuration  whenever  all  the  control  points  are 
at  their  goal  locations.)  The  goal  configuration  may  not  be  uniquely  defined.  For  instance, 
if  A  is  a  2D  object  that  can  both  translate  and  rotate  in  the  plane  (3D  configuration  space), 
the  specification  of  the  goal  positions  of  two  control  points  uniquely  determines  the  goal 
configuration  of  the  robot.  If  A  is,  say,  a  10-DOF  manipulator  arm,  specifying  the  desired 
positions  of  some  points  in  the  end-effector  determines  a  goal  region  in  configuration  space. 
In  all  cases,  we  denote  by  Cgoai  the  subset  of  goal  configurations. 

It  is  important  that  a  path  planner  allows  to  specify  a  goal  configuration  by  specifying  the 
goal  positions  of  a  small  number  of  points  in  the  robot.  Indeed,  in  many  tasks,  the  goal 
configuration  is  incompletely  determined  by  the  task  constraints.  Arbitrarily  selecting  one 
configuration  among  the  various  possible  ones  may  result  in  a  more  difficult,  even  impossible, 
path  planning  problem.  Furthermore,  if  the  robot  has  many  DOF's,  specifying  a  unique 
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goal  configuration  is  a  difficult  task  in  itself,  since  it  requires  a  collision-free  placement  of 
the  various  bodies  of  the  robot  to  be  found. 


In  order  to  precisely  define  the  C-potential,  we  must  now  specify  the  arbitration  function 
G.  In  most  of  the  previous  systems  using  the  artificial  potential  field  approach,  G  has  been 
chosen  as  a  positive  linear  combination  of  the  W-potentials  [Khatib,  1986]: 

i—k 

G(yi  j  •  •  • ,  yk)  =  ^  ^ i yi ♦ 

t=i 

This  simple  choice  seems  natural  because  it  does  not  favor  one  control  point  over  the  others. 
However,  precisely  for  that  reason,  conflicts  among  the  points  tend  to  be  frequent,  producing 
numerous  spurious  local  minima. 

The  choice  of  the  function  G  is  specially  important  because  it  highly  influences  the  number 
and  the  depths  of  the  local  minima  of  U.  With  W-potentials  free  of  local  minima,  the 
workspace  concavities  do  not  directly  create  local  minima  of  U.  It  is  the  concurrent  attrac¬ 
tion  of  the  different  control  points  toward  their  respective  goal  positions,  which  creates  these 
local  minima.  This  results  from  the  fact  that  these  points  do  not  move  independently  from 
each  other.  As  discussed  above,  the  function  G  determines  how  the  competition  between 
the  different  points  is  to  be  regulated. 

The  choice  of  G  which  seems  to  minimize  the  number  of  local  minima  is: 


i—k 


G(yi,...,yk)  =  rriiny;. 

2=1 


Indeed,  this  arbitration  function  favors  the  attraction  of  the  point  which  is  already  in  the 
best  position  to  reach  its  goal.  However,  when  one  point  has  reached  its  goal  position, 
the  potential  field  is  identically  zero,  and  the  other  points  are  not  attracted  toward  their 
goal  positions.  A  way  to  avoid  this  shortcoming  is  to  add  another  term  to  the  arbitration 
function: 


i—k 


i—k 


G{y\ ,  •  •  • ,  J/jfc)  =  min  yi  +  e  max  yi 
1  =  1  2  =  1 


(1) 


where  e  is  a  small  real  number.  In  our  experiments,  we  used  e  =  0.1.  However,  the  best 
value  of  e  may  depend  on  the  robot. 


Another  choice  for  G  is: 

i=k 

G(y\,.  ..,yk)  =  max  yi .  (2) 

This  arbitration  function  favors  the  attraction  of  the  control  point  that  is  the  furthest 
away  from  its  goal  (along  the  path  determined  by  the  W-potential).  It  tends  to  increase 
the  number  of  competitions  among  the  control  points  and,  therefore,  the  number  of  local 
minima.  However,  it  can  be  a  good  choice  for  robots  with  many  DOF’s.  As  a  matter 
of  fact,  the  number  of  local  minima  is  not  the  only  measure  for  the  quality  of  the  C- 
potential.  Another  critical  factor  is  the  depths  of  these  minima.  Indeed,  in  some  cases, 
it  may  be  preferable  to  have  several  small  local  minima,  rather  than  a  single,  very  deep 
one.  This  is  specially  true  when  the  robot  has  many  DOF’s,  since  the  number  of  discretized 
configurations  contained  in  a  local  minimum  well  of  given  depth  increases  exponentially  with 
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the  dimension  of  the  configuration  space.  Various  experiments  with  the  above  arbitration 
function  indicate  that  in  general  it  increases  the  number  of  local  minima,  but  reduces  their 
volumes. 

Obviously,  many  other  arbitration  functions  can  be  imagined.  In  some  of  our  experiments, 
we  used  other  C-potentials,  as  indicated  further  in  the  paper. 

We  now  describe  four  path  planning  techniques  incorporating  the  general  planning  approach 
presented  in  the  introduction.  All  these  techniques  iteratively  consider  each  level  of  reso¬ 
lution  in  the  workspace  and  configuration  space  pyramids,  from  the  coarsest  to  the  finest. 
They  terminate  with  success  as  soon  as  a  path  has  been  generated.  They  return  failure  if, 
after  having  considered  the  finest  bitmaps,  they  still  have  not  generated  a  path.  Below,  we 
only  describe  the  algorithms  executed  at  each  level  of  resolution.  We  denote  by  QW  and 
QC  the  current  workspace  and  configuration  space  grids. 

4  Best-First  Motion  Technique 

4.1  Description 

We  start  with  a  very  simple  path  planning  algorithm7.  This  algorithm  essentially  performs 
a  best-first  search  [Nilsson,  1980]  of  the  corresponding  configuration  space  discretized  grid 
QC  using  the  C-potential  U  as  the  cost  function  to  minimize.  This  means  that  the  possible 
successors  of  a  configuration  are  all  its  ^-neighbors,  for  some  k  6  [l,n],  in  the  grid  and 
these  neighbors  are  selected  by  decreasing  order  of  the  C-potential.  A  configuration  has  2 n 
1-neighbors,  2n2  2-neighborhoods,...,  and  3n  —  1  n-neighbors.  In  our  implementation  of  the 
best-first  motion  technique,  we  used  k  =  n.  The  size  of  the  neighborhood  thus  increases 
exponentially  with  n,  but  as  we  will  see  below,  the  planing  technique  is  intrinsically  limited 
to  planning  problems  involving  few  DOF’s  (typically,  n  <  4),  so  that  the  size  of  the  maximal 
neighborhood  remains  reasonable. 

As  long  as  the  algorithm  does  not  reach  a  local  minimum  of  the  C-potential,  the  search 
reduces  to  following  an  approximation  of  the  negated  gradient  of  the  C-potential  (fastest 
descent  procedure).  When  a  local  minimum  is  reached,  the  same  algorithm  naturally  fills 
up  the  local  minimum  well  until  a  saddle  point  is  reached.  Then,  it  proceeds  again  along 
the  negated  gradient  of  U.  It  stops  when  a  goal  configuration  (U  =  0)  is  attained.  Hence, 
the  algorithm  basically  searches  a  graph  connecting  local  minima  among  them,  but  this 
graph  is  not  explictly  represented.  Between  two  minima,  it  simply  executes  a  deterministic 
gradient  motion.  At  each  local  minimum  it  spends  so  time  filling  it  up.  Two  local  minima 
are  adjacent  in  the  local  minima  graph  if  they  communicates  through  a  saddle  point.  From 
a  local  minimum,  the  algorithm  always  goes  to  the  adjacent  local  minimum  attainable 
through  the  lowest  saddle  point.  “Backtracking”  occurs  in  the  local  minima  graph  when 
two  adjacent  minima  get  filled  up;  then,  they  are  implicitly  merged  into  a  single  minimum, 
their  common  filling  level  being  the  new  height  of  the  minimum. 

The  algorithm  is  resolution-complete,  i.e,  it  is  guaranteed  to  reach  the  goal  in  a  finite 

amount  of  time  whenever  a  solution  path  exists  (at  the  maximal  resolution),  or  return 

7This  algorithm  is  described  in  more  detail  in  [Barraquand  and  Lat.ombe,  1989]. 
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failure  when  there  is  no  solution.  In  most  cases  (virtually  all  reasonable  cases),  a  very 
small  subset  of  the  grid  QC  has  to  be  explored  before  the  algorithm  terminates.  It  is  easy, 
however,  to  create  a  problem  admitting  no  solution  path,  which  requires  most  of  QCfree  to 
be  explored  before  the  planner  gives  up. 

At  this  stage,  there  is  an  important  aspect  of  the  algorithm  that  must  be  made  precise. 
Since  the  C-potential  only  depends  on  the  distance  of  a  few  control  points  to  the  obstacles, 
the  best-first  search  does  not  guarantee  that  collisions  of  the  robot  with  the  obstacles  will 
be  completely  avoided.  Therefore,  whenever  the  planner  considers  a  new  configuration  in 
QC,  it  should  check  that  it  lies  in  the  free  space.  Because  the  planner  does  not  represent 
the  C-obstacles  explicitly,  the  verification  is  done  in  the  workspace  as  explained  below. 

As  the  workspace  representation  is  distributed,  it  seems  natural  to  represent  the  robot  A 
in  the  same  way,  i.e.  by  a  bitmap.  If  we  were  using  a  massively  parallel  computer,  this 
representation  would  be  the  simplest  and  the  most  efficient  to  test  collisions.  For  example, 
we  could  have  one  processor  per  point  in  the  bitmap  representation  of  A;  this  processor 
would  compute  the  position  of  the  point  in  the  workspace  at  the  current  configuration  of 
A  and  determine  whether  it  is  contained  in  the  obstacle  region,  or  not.  However,  as  we 
implemented  our  planner  on  a  sequential  computer,  we  chose  to  represent  the  boundary 
of  A  algebraically..  The  implemented  collision  checking  technique  is  a  classical  “divide- 
and-conquer”  algorithm,  which  consists  of  iteratively  adjusting  the  number  of  points  in  the 
robot’s  boundary  which  are  necessary  to  check  collisions  reliably.  To  illustrate  this  idea, 
let  us  consider  a  simple  robot  modelled  as  a  single  straight  line  segment  of  length  L.  The 
function  d\  has  already  been  computed  over  Wempty  Instead  of  checking  the  distance  di  for 
each  discretized  point  on  the  line  segment  modelling  A  and  then  calculating  the  minimum 
of  all  these  distances,  we  first  compute  the  distances  di(begin(q ))  and  di(end(q))  of  the 
two  extremity  points  begin  and  end  of  the  segment.  If  min{di(i>eym(q)),  di(end(q))}  >  L, 
then  we  are  certain  that  the  robot  does  not  collide  any  obstacle.  Otherwise,  we  cut  the 
segment  in  two  half-length  segments,  and  we  recursively  apply  the  procedure  to  the  two 
new  segments.  If  the  robot  boundary  is  modelled  by  a  polygon,  the  same  procedure  can  be 
applied  to  each  edge  of  the  polygon.  The  procedure  can  also  be  generalized  to  higher-order 
shapes  of  the  robot  boundary  such  as  circular  and  elliptical  arcs.  When  the  robot  is  far 
from  the  obstacles,  very  few  points  on  its  boundary  need  to  be  checked.  When  the  robot 
gets  closer  to  an  obstacle,  an  increasing  number  of  points  (up  to  the  current  resolution  of 
the  workspace  bitmap)  in  the  boundary  segment  that  is  close  to  make  contact  have  to  be 
checked . 

The  search  algorithm  only  considers  the  neighbors  of  a  configuration  that  are  collision- 
free.  Hence,  there  may  be  two  types  of  local  minima:  the  natural  minima  of  U  (where  the 
gradient  is  zero)  and  the  minima  located  at  the  boundary  of  QCfree  (where  the  gradient 
is  not  zero  in  general).  When  an  obstacle  is  hit,  the  last  configuration  before  the  collision 
(i.e.,  in  the  boundary  of  QCfree )  is  taken  as  the  local  minimum.  Both  types  of  minimum 
are  escaped  in  the  same  fashion. 
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Figure  9:  Path  Generated  for  a  3-DOF  Mobile  Robot 
4.2  Experimental  Results 

We  experimented  the  path  planner  on  a  “mobile  robot”  with  two  DOF’s  of  translation  and 
one  DOF  of  rotation,  namely  a  long  rectangular  bar  in  a  2D  workspace.  Figure  9  shows 
an  example  of  path  generated  by  the  planner.  This  path  demonstrates  the  ability  of  the 
planning  technique  to  produce  complex  maneuvers.  The  running  time  to  produce  the  path 
was  1  second8.  This  is  three  orders  of  magnitude  faster  than  the  running  times  reported 
in  [Brooks  and  Lozano-Perez,  1983]  for  similar  (though  apparently  simpler)  path  planning 
problems.  One  order  of  magnitude  is  due  to  the  faster  computer  we  used.  The  other  two 
are  actually  a  product  of  our  algorithm. 

Figure  10  shows  another  example  of  the  planning  abilities  of  the  algorithm.  For  the  same 
robot,  a  path  was  generated  in  less  than  5  seconds  within  a  5122  bitmap  representing  a 
workspace  cluttered  by  more  than  70  complex-shaped,  randomly  generated  obstacles.  This 
example  would  probably  be  very  difficult  to  run  with  a  planner  using  a  semi-algebraic 
representation  of  the  workspace,  e.g.  an  exact  cell  decomposition  planning  technique.  It 
demonstrates  the  power  of  the  “distributed”  representations  used  in  our  planner  relatively 
to  the  “centralized”  algebraic  representations. 

In  both  examples,  the  C-potential  was  computed  by  considering  two  control  points  located 
at  the  midpoints  of  the  two  small  edges  at  the  extremities  of  the  bar,  using  the  improved  W- 
potential  described  in  Subsection  3-2.2.  The  C-potential  was  computed  using  the  arbitration 
function  defined  by  formula  (1).  At  first  sight,  the  experimental  efficiency  of  the  best-first 

8 Since  there  may  be  many  simple  variants  of  the  techniques  presented  in  this  paper,  only  the  order  of 
magnitude  of  the  given  execution  times  is  actually  pertinent. 
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Figure  10:  Path  Generated  Among  Randomly  Distributed  Obstacles 

motion  technique  is  surprising.  The  fundamental  reason  of  this  efficiency  is  the  following. 
The  potential  fields  computed  in  the  workspace  are  designed  to  be  “perfect”  potential  fields 
for  a  point  robot,  in  the  sense  that  they  have  no  other  minimum  than  the  goal.  Therefore, 
complex  shaped  obstacles  with  large  concavities  do  not  directly  create  local  minima  of 
the  C-potential.  Local  minima  occur  only  where  the  workspace  is  so  cluttered  that  the 
solution  path  has  to  come  very  close  to  the  obstacles.  But  in  that  case,  the  local  minima 
are  close  to  the  boundary  of  the  free  space.  Hence,  they  have  small  domains  of  attraction, 
i.e.  contain  a  small  number  of  discrete  configurations.  The  algorithm  fills  the  wells  very 
quickly.  Furthermore,  because  it  guides  the  robot  along  a  path  where  maneuvering  space  is 
maximized,  the  improved  W-potential  tends  to  produce  fewer  local  minima.  The  algorithm 
also  works  in  a  satisfactory  fashion,  but  is  in  general  slower,  when  the  simple  W-potential 
is  used  in  place  of  the  improved  one. 

However,  in  practice,  the  above  planning  technique  is  only  applicable  to  robots  with  a  small 
number  n  of  DOF’s  -  typically,  n  <  4.  Indeed,  the  number  of  discrete  configurations  in  a 
local  minimum  well  increases  exponentially  with  the  number  of  DOF’s.  Filling  up  a  well 
would  be  too  time  consuming  for  robots  with  many  DOF’s. 

5  Random  Motion  Technique 

5.1  Description 

The  path  planning  algorithm9  presented  in  this  section  essentially  differs  from  the  previous 
one  in  the  way  it  escapes  local  minima.  Rather  than  filling  up  each  encountered  local  min- 

9This  algorithm  is  described  in  more  detail  in  [Barraquand  and  Latombe,  1989]. 
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imum,  it  applies  a  Monte-Carlo  procedure,  which  consists  of  generating  Brownian  motions 
until  the  minimum  is  escaped. 

Starting  from  the  initial  configuration,  the  algorithm  first  searches  the  current  grid  QC  in  a 
best-first  fashion.  Thus,  it  follows  the  negated  gradient  of  the  C-potential,  until  it  reaches 
a  local  minimum  (call  it  q/oc).  We  call  such  a  motion  a  gradient  motion.  If  q \oc  E  Cg0ah 
the  planner  returns  the  solution  path  generated.  Otherwise,  it  attempts  to  escape  the  local 
minimum  by  generating  several  random  motions  from  q/oc.  These  random  motions  are 
Brownian  motions  described  in  more  detail  in  the  next  subsection. 

The  generation  of  gradient  motions  and  the  recognition  of  local  minima  are  done  using  the 
n-neighborhood  in  QC.  As  long  as  n  <  4,  this  raises  no  difficulty.  However,  when  n  becomes 
too  big  (and  we  ran  the  random  motion  technique  for  robots  with  many  DOF’s),  the  size  of 
the  n-neighborhood  becomes  too  large.  In  this  case,  at  each  step  of  the  gradient  motion,  only 
a  few  points  in  the  n-neighborhood  of  the  current  position  x  are  considered.  This  is  done  by 
exploring  the  n-neighborhood  of  x  iteratively  in  a  random  order  and  limiting  the  number 
of  iterations.  At  each  iteration,  a  point  x'  is  randomly  chosen  in  the  n-neighborhood  of  x 
(using  a  uniform  distribution  law).  If  U(x')  <  U(x),  x'  is  taken  as  the  successor  of  x  along 
the  path  of  the  gradient  motion  (hence,  the  motion  may  only  follow  a  rough  approximation 
of  the  negated  gradient  flow).  The  number  of  iterations  is  limited  to  a  few  tens  to  a  few 
hundreds  (depending  on  the  value  of  n).  If  none  of  them  generates  a  successor  of  x,  x  is 
treated  as  a  local  minimum.  An  alternative  to  this  technique  would  be  to  use  a  smaller 
neighborhood,  for  instance  the  linear  1-neighborhood.  In  fact,  we  tried  several  alternatives, 
and  the  technique  described  above  gave  the  best  results.  In  particular,  due  to  the  crude 
discretization  that  it  entails,  a  small  neighborhood  often  resulted  in  the  detection  of  many 
fictitious  local  minima. 

At  the  terminal  configuration  of  every  random  motion,  the  algorithm  executes  a  gradient 
motion  until  it  reaches  an  hopefully  new  local  minimum.  From  each  local  minimum,  if  none 
of  them  is  in  the  goal  region,  it  performs  another  set  of  random  motions.  Etc...  The  graph  of 
the  local  minima  is  thus  incrementally  built,  the  path  joining  two  “adjacent”  local  minima 
being  the  concatenation  of  a  random  motion  and  a  gradient  motion.  Using  the  values  of  the 
C-potential  at  the  local  minima,  a  best-first  search  of  this  graph  is  performed  until  the  goal 
configuration  is  reached  or  all  the  local  minima  in  the  graph  have  been  expanded  (i.e.,  a  series 
of  random  motions  have  been  executed  from  them).  When  a  solution  path  is  generated, 
it  is  smoothed  using  a  classical  variational  calculus  technique.  An  interesting  property  of 
this  algorithm  is  that  all  the  random  motions  starting  from  a  given  local  minimum  can  be 
performed  concurrently  on  a  parallel  machine,  since  there  is  no  need  for  communication 
between  the  different  processing  units. 

Many  variants  of  the  best-first  search  scheme  are  possible,  and  we  experimented  with  several 
ones.  The  search  technique  that  gave  the  best  results  consists  of  iteratively  executing 
random  motions  starting  at  the  current  local  minimum  q \oc  and,  after  each  one,  performing 
a  gradient  motion.  If  the  attained  local  minimum  has  a  lower  C-potential  than  q/oc,  the 
iteration  is  stopped  and  the  search  proceeds  from  this  new  local  minimum,  by  executing 
new  random  motions.  The  number  of  random  motions  starting  at  each  local  minimum  is 
arbitrarily  bounded.  If  no  random  motion  from  q/oc  followed  by  a  gradient  motion  leads 
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the  planner  to  a  better  minimum,  q ioc  is  considered  as  a  dead-end.  The  algorithm  then 
backtracks  to  a  point  in  the  current  path  connecting  the  initial  configuration  to  q \oc.  This 
point  is  selected  randomly,  using  a  uniform  distribution  law  over  the  set  of  points  contained 
in  the  current  path.  The  search  for  a  path  resumes  at  that  point  by  executing  a  gradient 
motion.  Since  the  backtracking  point  may  belong  to  the  subpath  generated  by  a  random 
motion,  the  gradient  motion  may  lead  to  a  local  minimum  not  encountered  yet.  This  search 
technique  has  the  advantage  of  not  requiring  an  explicit  representation  of  the  local  minima 
graph  to  be  maintained. 

As  the  algorithm  uses  a  random  procedure  to  build  the  graph  of  the  local  minima,  it  is 
not  guaranteed  to  find  a  path  whenever  one  exists.  In  other  words,  the  algorithm  is  not 
complete.  However  the  properties  of  Brownian  motions  make  it  possible  to  prove  that  when 
the  number  of  Brownian  motions  executed  from  every  local  minimum  is  unbounded  (the 
computation  time  may  tend  towards  infinity),  the  probability  to  reach  the  goal  converges 
towards  1  (see  [Barraquand  and  Latombe,  1989]).  Hence,  we  say  that  the  algorithm  is 
probabilistically  resolution-complete.  However,  this  convergence-in-distribution  prop¬ 
erty,  which  is  well-known  for  the  so-called  “simulated  annealing”  algorithms10  (see  [Geman 
and  Geman,  1986]),  is  a  very  weak  one.  Indeed,  the  totally  uninformed  algorithm,  which 
executes  a  Brownian  motion  from  q ina  and  terminates  when  it  enters  a  small  neighborhood 
of  the  goal  configuration,  is  also  probabilistically  resolution-complete! 

A  weakness  of  the  algorithm  with  the  search  technique  described  above  is  that  it  may  be 
unable  to  recognize  that  a  planning  problem  admits  no  solution  path.  However,  since  the 
algorithm  turns  out  to  be  relatively  fast  in  most  practical  cases  where  there  exist  solution 
paths  (see  Subsection  5.3),  one  way  to  deal  with  this  drawback  is  to  arbitrarily  limit  the 
computing  time  allocated  to  the  planner.  By  simplifying  a  bit,  we  can  say  that  we  traded 
decidability  and  slow  computation  for  semi-decidability  and  fast  computation. 

5,2  Brownian  Motions 

In  theory,  a  Brownian  motion  is  a  continuous  stochastic  process.  In  our  planner,  it  is 
implemented  as  a  discrete  random  walk.  It  consists  of  discretizing  time  into  units  and 
executing  a  motion  step,  whose  projection  along  each  axis  <?*,  i  =  1,  ...,n,  is  randomly  A; 
or  —At‘,  at  each  time  unit.  The  motion  lasts  a  certain  number  of  time  units,  say  t.  At  each 
time  unit,  the  motion  step  is  independent  of  the  previously  executed  steps.  This  random 
walk  is  known  to  converge  almost  surely  toward  a  Brownian  motion  when  every  A t*  tends 
toward  0  [Papoulis,  1965]. 

Without  lack  of  generality,  let  us  take  the  local  minimum,  q/oc,  as  the  origin.  The  con¬ 
figuration  attained  by  a  Brownian  motion  of  duration  t  is  a  random  variable  Q (t)  = 
(Qi(t), ...,  Qn(t)),  with  the  following  properties  [Papoulis,  1965]: 

-  The  density  Pi{qi)  of  Qi(t)  is: 

1  2 

i0For  a  comparison  between  simulated  annealing  algorithms  and  our  planning  algorithm,  see  [Barraquand 
and  Latombe,  1989]. 
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-  The  standard  deviation  D{  of  Qi(t)  is: 

Di  =  A  iVt. 

The  Brownian  motion  is  well-defined  as  long  as  it  does  not  encounter  any  obstacle  in  con¬ 
figuration  space.  When  it  hits  the  boundary  of  a  C-obstacle,  the  Brownian  motion  has 
to  be  adapted  so  that  it  stays  in  the  free  space.  The  classical  adaptation  of  a  Brownian 
motion  when  the  space  is  bounded  consists  of  reflecting  the  motion  that  would  take  place 
if  there  were  no  boundary,  symmetrically  to  the  tangent  hyperplane  of  the  boundary  at 
the  collision  configuration.  The  mathematical  consistency  of  this  adaptation  is  discussed 
in  detail  in  [Anderson  and  Orey,1976].  Our  planner,  which  does  not  construct  an  explicit 
representation  of  the  C-obstacles,  does  not  know  the  orientation  of  the  tangent  hyperplane 
at  the  collision  configuration.  Hence,  whenever  a  motion  step  leads  to  a  collision,  instead  of 
bouncing  symmetrically  on  the  hyperplane  tangent  to  the  C-obstacle,  the  planner  guesses 
another  random  step  and  substitutes  it  for  the  previous  one.  The  collision  is  detected  in 
the  workspace  using  the  divide-and-conquer  technique  described  in  Subsection  4.1. 


The  duration  t  of  a  Brownian  motion  still  has  to  be  chosen.  If  it  is  too  short,  the  motion 
has  little  chance  to  escape  the  local  minimum.  If  it  is  too  long,  the  planner  may  waste  time 
and  loose  the  opportunity  of  using  the  C-potential  gradient  information  when  it  becomes 
useful  again.  We  can  define  the  attraction  radius  Am{<lloc )  of  any  local  minimum  q/oc  of 
U  along  each  axis  qi  as  the  distance  along  qi  between  q/oc  and  the  nearest  saddle  point  of 
U  in  that  direction.  In  order  to  escape  the  local  minimum  q/oc,  the  minimum  distance  that 
the  robot  has  to  travel  in  each  direction  qi  from  q/oc  is  precisely  A#z*(q/oc).  If  we  were  able 
to  estimate  the  statistics  of  Arx ,  the  property  D{  =  A t\/T  would  give  us  a  clue  for  escaping 
the  minimum.  The  duration  of  the  Brownian  motion  would  be  the  function  t(q/oc)  defined 
by: 


t(q ioc)  ~  max 

*€[l,n] 


/  ARi(qioc) 

V  A, 


) 


2 


(3) 


But,  as  we  make  no  assumption  on  the  obstacle  distribution,  we  cannot  infer  any  strong 
statistical  property  about  U  and  Anr  However,  in  general,  we  may  consider  that  the 
distance  Am  for  each  parameter  qi  does  not  exceed  the  distance  that  would  provoke  a 
motion  of  the  robot  longer  than  the  workspace  diameter  itself.  Normalizing  this  diameter 
to  1,  this  assumption  yields  the  following  estimate  of  Am  for  any  local  minimum  q/oc: 


Am  &  l/JslUp  (4) 

This  estimate,  combined  with  the  fact  that  Ar{  is  strictly  positive,  leads  us  to  treat  Arx  as 
a  random  variable  with  a  truncated  Laplace  distribution  of  density: 


p{ARi)  =  Js'upexp(-Jstupylfit) 

(The  truncated  Laplace  distribution  is  the  less  informed  distribution  -  i.e.  the  one  which 
maximizes  entropy  -  for  a  positive  random  variable  of  given  expected  value.) 


We  have  A;  ~  <5/^sup  (see  Subsection  2.3),  where  6  is  the  distance  between  two  adjacent 
points  in  QW .  Combining  this  formula  with  (3)  and  (4),  we  obtain: 
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Figure  11:  Path  Generated  for  a  3-DOF  Robot 


One  could  take  t  equal  to  the  above  value.  However,  this  choice  would  implicitly  assume 
that  all  the  attraction  radii  are  the  same,  which  is  not  the  case.  Using  the  above  distribution 
of  Ar{,  we  choose  t  as  the  value  of  a  random  variable  T  with  the  following  density: 


P(t)  = 


6_ 

2  sft 


exp(— 8y/t). 


One  can  verify  that  the  expected  value  of  T  is  indeed  l/<52. 

In  fact,  a  value  of  T  gives  a  maximal  duration  of  the  Brownian  motion.  At  each  time  unit 
during  the  motion,  the  planner  checks  the  C-potential  at  the  current  configuration  against 
U(q /oc).  If  it  is  smaller,  the  planner  terminates  the  Brownian  motion. 


5.3  Experimental  Results 

We  have  tested  the  random  motion  technique  with  several  different  robot  stuctures.  We 
describe  now  some  of  the  obtained  results.  Since  the  planning  technique  contains  several 
random  components,  the  running  time  for  the  same  planning  problem  is  not  constant.  The 
times  given  below  are  only  typical  times.  It  is  not  unusual  that  two  running  times  for  the 
same  example  differ  in  a  ratio  of  5.  The  execution  times  would  be  both  smaller  and  much 
more  stable  on  a  parallel  architecture  allowing  the  concurrent  execution  of  several  random 
motions. 

First,  we  applied  it  to  the  3-DOF  robot  example  shown  in  Figure  9  using  the  same  C- 
potential  as  with  the  best-first  motion  algorithm.  We  got  the  path  shown  in  Figure  11  (after 
smoothing).  The  overall  computation  time  was  10  seconds,  which  is  about  ten  times  slower 
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Figure  12:  Structure  of  the  8-DOF  Manipulator 

than  with  the  pure  best-first  search.  However,  we  think  that  a  parallel  implementation  of 
the  random  motion  technique  will  considerably  reduce  this  computation  time. 

We  experimented  with  the  planner  on  a  8-DOF  serial  manipulator  with  only  revolute  joints 
(see  Figure  12),  Figure  13  illustrates  a  path  generated  by  the  planner.  The  goal  region  is 
defined  by  the  position  of  the  end-point  p  of  the  last  link  of  the  robot.  The  C-potential  U 
was  computed  as  the  improved  W-potential  at  p ,  i.e.  U p(X(p,  q)).  The  path  was  generated 
in  2  minutes,  covering  all  the  computation  performed  between  the  input  of  the  robot  and 
workspace  models  and  the  display  of  the  generated  path. 

We  also  experimented  with  the  planner  using  a  10-DOF  non-serial  manipulator  arm  with 
3  prismatic  and  7  revolute  joints,  as  depicted  in  Figure  14.  Figures  15  illustrates  a  path 
generated  by  the  planner.  The  C-potential  was  computed  by  considering  two  control  points 
located  at  the  end-points  of  the  two  kinematic  chains,  using  the  improved  W-potential  field 
and  the  arbitration  function  defined  in  formula  (2).  The  computation  time  was  3  minutes. 
(In  this  example,  the  number  of  configurations  in  the  grid  QC  in  which  the  path  was  found 
is  of  the  order  of  lOO10  =  lO20.) 

Finally,  we  have  run  the  path  planning  algorithm  with  a  31-DOF  serial  manipulator  carrying 
a  bar  in  a  3D  workspace.  The  manipulator  (see  Figure  16)  consists  of  a  sequence  of  10 
identical  modules,  each  with  three  DOF’s  (one  translation  and  two  rotations).  The  last 
module  provides  a  third  rotation.  A  path  generated  by  the  planner  is  illustrated  in  Figure 
17.  The  C-potential  was  computed  by  considering  two  control  points  located  at  the  end¬ 
points  of  the  bar,  using  the  improved  W-potential  field  and  the  arbitration  function  of 
formula  (2).  The  computation  time  was  15  minutes.  The  level  of  resolution  of  the  workspace 
representation  required  to  find  this  path  was  I283. 
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Figure  14:  Structure  of  the  10-DOF  Manipulator 
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Figure  16:  Structure  of  the  31-DOF  Manipulator 
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ure  17:  Path  Generated  for  the  31-DOF  Manipulator 


In  all  the  above  multi-link  robot  examples,  we  simulated  mechanical  stops  by  limiting  the 
range  of  displacement  of  every  joint.  Except  with  the  31-DOF  arm,  we  also  prohibited 
collisions  among  the  links,  by  checking  collisions  between  any  two  links  using  the  technique 
of  Subsection  4.1. 

Again  the  efficiency  of  the  algorithm  may  be  surprising.  We  think  that  it  results  from  the 
fact  that  a  typical  path  planning  problem  has  many  solutions,  so  that  a  globally  random 
process  can  find  one,  providing  it  is  sufficiently  informed  most  of  the  time.  This  efficiency 
result  is  not  new.  Monte-Carlo  techniques  have  already  been  helpful  for  solving  NP-complete 
problems.  For  example,  Cerny  [Cerny,  1985]  described  an  algorithn  that  generates  sub- 
optimal  solutions  to  the  traveling  salesman  problem  for  more  than  10,000  towns.  Kirkpatrick 
et  al.  [Kirkpatrick,  Gelatt  and  Vecchi,  1983]  proposed  an  algorithm  for  the  placement  and 
the  routing  of  VLSI  chips,  which  is  better  than  human  experts.  In  both  problems,  the  very 
large  search  space  is  associated  with  a  large  number  of  “good”  sub-optimal  solutions. 

6  Valley-Guided  Motion  Technique 

6.1  Description 

The  path  planning  algorithm11  described  in  this  section  differs  significantly  from  the  previ¬ 
ous  two.  It  consists  of  searching  the  “valleys”  of  a  C-potential  U  in  Cfree.  The  set  of  valley 
points  of  U  is  called  the  valley  roadmap  and  is  denoted  by  V. 

The  algorithm  is  similar  in  structure  to  a  retraction  algorithm  [O’Dunlaing,  Sharir  and  Yap, 
1983]: 

1.  Given  an  initial  and  a  goal  configurations,  q inn  and  q goai,  generate  two  new 
configurations,  q,  and  q^,  local  minima  of  U,  by  following  the  negated  gradient 
of  U  from  qinn  and  qgoaL  respectively. 

2.  Search  the  valley  roadmap  V  for  a  path  connecting  q*  to  qg. 

3.  If  step  2  terminates  successfully,  return  the  path  obtained  by  concatenating  the 
three  paths  joining  respectively  q to  q ;  (gradient  motion),  q t  to  q^  (valley- 
guided  motion),  and  q^  to  q goal  (gradient  motion).  Otherwise,  return  failure. 

This  algorithm  does  not  require  U  to  be  minimum  at  the  goal  configuration.  However,  since 
<lgoal  is  projected  in  V  by  a  gradient  motion,  qgoai  has  to  be  uniquely  defined,  which  was  not 
the  case  with  the  previous  two  algorithms12.  Ideally,  U  should  be  such  that  V  is  a  ID  subset 
of  the  free  space  C/ree.  In  addition,  U  should  “represent  the  topology”  of  C/ree,  i.e.  q *  and 
q5  should  be  connected  in  V  iff  q init  and  q goai  are  connected  in  Cfree.  In  [Barraquand, 
Langlois  and  Latombe,  1989],  we  discuss  the  conditions  under  which  a  C-potential  satisfies 

11  This  algorithm  is  described  in  more  detail  in  [Barraquand,  Langlois  and  Latombe,  1989], 

12Actually,  one  could  extend  the  above  algorithm  to  handle  the  case  where  a  goal  configuration  is  any 
configuration  in  a  subset  Cgoai  of  the  configuration  space.  But  this  would  require  each  configuration  in  the 
discretized  boundary  of  Cgoai  to  be  connected  by  a  gradient  motion  to  V.  Except  for  robots  with  few  DOF’s, 
this  is  likely  to  be  impractical. 
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these  two  properties.  Unfortunately,  these  conditions  are  quite  involved  and  difficult  to 
verify. 

At  step  2,  V  is  searched  in  a  depth-first  manner.  A  decision  is  made  at  every  crossroad 
using  a  simple  heuristic  function  defined  as  another  C-potential  Ufeeiir.  This  C-potential  is 
constructed  as  described  in  Section  3,  i.e.  by  combining  local-minima-free  W-potentials.  In 
our  experiments,  this  heuristic  potential  dramatically  reduced  the  computation  times. 

It  now  remains  to  specify  how  V  is  constructed.  Unfortunately,  as  noticed  in  [Barraquand, 
Langlois  and  Latombe,  1989],  there  has  not  been  much  research  on  the  concept  of  valley 
reported  in  the  literature  so  far.  One  way  to  get  a  simple,  but  ad  hoc,  definition  of  val¬ 
ley  points  is  to  discretize  the  set  of  possible  directions  of  the  valleys.  Namely,  we  force 
these  directions  to  be  the  various  qi  coordinate  axes.  Thus,  we  only  have  n  possible  valley 
directions.  We  define  a  configuration  q  to  be  a  valley  point  along  the  qi  coordinate  axis 
(*’  €  [l,n])  when  all  the  values  of  the  C-potential  U  at  the  (2n  -  2)  1-neighbors  in  the 
hyperplane  orthogonal  to  this  direction  are  greater  than  U(q).  According  to  this  simple 
definition,  every  local  minimum  of  U  is  also  a  valley  point  of  U  (which  is  consistent  with 
our  intuitive  notion  of  valley),  so  that  the  valley  roadmap  V  can  be  regarded  as  a  graph 
connecting  the  local  minima  of  U.  This  graph  has  to  be  searched  for  a  path  between  two 
minima,  q,  and  qff. 

In  order  to  check  whether  a  point  q  is  a  valley  point  or  not,  we  thus  use  the  simple  following 
algorithm: 

1.  Compute  U(q). 

2.  Compute  the  2n  values  of  U  at  the  1-neighbors  of  q. 

3.  For  each  possible  valley  direction  i  £  [l,n]  do: 

Compare  U(q)  to  the  2n  —  2  values  of  U  at  the  1-neighbors  in  the  hyperplane 
orthogonal  to  the  qi  axis.  If  U(q)  is  smaller  or  equal  to  these  2n  -  2  values,  q  is 
a  valley  point. 

The  time  complexity  of  this  algorithm  is  0(n2). 

As  non-degenerated  valleys  are  ID,  the  safest  neighborhood  to  track  them  is  the  n- 
neighborhood.  Unfortunately,  the  cardinal  of  this  neighborhood  is  3n  -  1,  yielding  an 
exponential  time  tracking  algorithm.  For  instance,  in  a  10D  configuration  space,  using  the 
10-neighborhood  entails  performing  about  60,000  valley  checkings  at  every  increment  of  the 
tracking  process.  Instead,  the  implemented  planner  uses  the  2-neighborhood,  whose  size  is 
quadratic  in  n,  allowing  each  tracking  increment  to  be  performed  in  0(n4)  time.  This  com¬ 
promise  between  time  and  completeness  appeared  reasonable  in  most  of  the  experiments  we 
made.  Complex  paths  have  been  generated,  meaning  that  valleys  are  trackable  using  the 
2-neighborhood.  Nevertheless,  the  choice  of  the  2-neighborhood  is  empirical. 

Some  valleys  may  be  dead-ends.  Valley  tracking  stops  when  the  tracking  algorithm  is  unable 
to  proceed  further.  It  also  stops  if  the  valley  is  running  into  an  obstacle  in  configuration 
space  (this  can  be  detected  by  using  the  collision  checking  technique  described  in  Subsection 
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Figure  18:  Path  Generated  for  a  3-DOF  Mobile  Robot 


4.1).  In  both  cases,  the  planning  algorithm  backtracks  to  a  crossroad  point  in  the  valley 
roadmap. 

6 . 2  Exp  erimentation 

We  experimented  with  the  valley-guided  motion  technique  on  a  variety  of  path  planning 
problems13.  The  best  results  were  obtained  with  the  C-potential  U  defined  as  follows 
[Barraquand,  Langlois  and  Latombe,  1989]: 

U(q)  =  <rlog£exp(UPi(X(Pi,q))/<T)) 

t  =  l 

where  cr  is  a  small  number  and  the  s  (i  =  l,...,s)  are  the  control  points14.  Each  W- 
potential  Up.  is  computed  using  the  simple  definition  of  Subsection  3.2.1. 

The  heuristic  potential  U heur  used  in  our  experiments  is  defined  as  the  sum  of  the  simple 
W-potentials  computed  at  a  few  control  points. 

We  applied  the  algorithm  to  3-DOF  mobile  robot  problems.  An  example  of  generated  path 
for  a  triangular  robot  is  shown  in  Figure  18.  This  example  was  run  using  15  control  points 

13  We  developed  the  valley-guided  motion  technique  before  the  three  other  techniques  presented  in  this 
paper.  It  has  been  implemented  on  a  SUN  III  workstation,  and  the  experiments  described  below  were  carried 
out  on  this  machine. 

14The  expression  defining  U(q)  is  a  discrete  approximation  of  a  log(J^  exp(Up(A(p,  q))/<r)dp),  which 
decreases  with  the  distance  separating  the  robot  and  the  obstacles  and  tends  toward  infinity  when  this  dis¬ 
tance  tends  toward  0.  The  expression  defining  U(q)  also  uniformly  converges  toward  the  simpler  expression 
sup-€[1  a]  Up(X(p, q))  when  <x  0.  The  latter  expression,  however,  is  highly  degenerated  for  manipulator 
arms  and  produces  flat  valleys  which  are  difficult  to  track.  Therefore,  it  cannot  be  used  reliably  by  the 
valley-guided  planning  technique. 
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Figure  19:  Path  Generated  for  the  10-DOF  Manipulator 

equally  distributed  in  the  triangle  boundary.  We  ran  the  algorithm  with  and  without  the 
heuristic  potential  U heur-  When  used,  XJheur  was  computed  by  considering  a  single  control 
point.  In  both  cases,  the  program  found  a  solution  path.  However,  the  heuristic  potential 
considerably  reduced  the  overall  computational  time.  In  the  example  of  Figure  18,  the 
running  time  was  reduced  from  2  minutes  to  less  than  30  seconds. 

We  also  applied  the  algorithm  to  the  10-DOF  robot  shown  in  Figure  14.  Figure  19  shows 
a  path  generated  by  the  planner.  In  this  example,  the  C-potential  U  was  computed  using 
70  control  points  equally  distributed  in  the  robot  and  the  heuristic  potential  U heuT  using 
two  control  points  located  at  the  end-points  of  the  two  kinematic  chains.  The  use  of  the 
heuristic  potential  was  necessary  to  obtain  a  path  in  a  reasonable  amount  of  time.  (In  a  10D 
configuration  space,  each  point  has  200  2-neighbors.  The  computation  of  the  next  point 
along  a  valley  takes  a  few  seconds.)  The  path  was  generated  in  less  than  1  hour. 

Notice  that  the  path  planning  problem  shown  in  Figure  19  is  significantly  simpler  than  the 
problem  shown  in  Figure  15.  There  is  more  space  for  the  robot  to  “maneuver”.  (In  Figure 
15,  the  hole  in  the  wall  is  narrower  and  the  square  object  is  closer  from  the  wall.)  The 
above  algorithm  failed  to  solve  the  problem  of  Figure  15.  The  search  was  stopped  after  a 
few  hours. 

The  above  planning  technique  is  slower  than  the  best-first  and  random  motion  techniques 
for  robots  with  few  DOF’s.  It  is  also  slower  and  definitely  less  reliable  than  the  random 
motion  technique  for  robots  with  many  DOF’s.  The  failures  of  the  algorithm  may  have 
several  causes:  the  search  of  the  valley  roadmap  may  require  exponential-time  computation, 
the  valley  roadmap  may  imperfectly  represent  the  connectivity  of  the  free  space,  it  may  be 
locally  degenerated  (i.e.,  may  not  be  ID),  the  2-neighborhood  may  be  insufficient  to  reliably 
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track  the  valleys,  etc...  We  think  that  these  causes  might  be  partly  remedied  by  investigating 
the  concept  of  valley  further  (see  [Barraquand,  Langlois  and  Latombe,  1989]). 

Notice  that  using  the  concept  of  valley  might  also  be  useful  to  improve  the  random  mo¬ 
tion  technique.  Indeed,  when  a  gradient  motion  reaches  a  local  minimum,  that  algorithm 
executes  a  series  of  non-informed  Brownian  motions.  Instead,  valley  directions  could  be 
extracted  around  the  encountered  minimum  and  the  probability  laws  governing  the  random 
motions  could  be  adapted  so  that  these  motions  tend  to  stay  closer  than  pure  Brownian 
motions  from  the  valley  directions. 

7  Constrained  Motion  Technique 

7.1  •  Description 

The  constrained  motion  technique  rests  on  a  combination  of  ideas  already  present  in  the 
previous  three  techniques.  On  the  one  hand,  like  the  best-first  and  random  motion  tech¬ 
niques,  it  makes  use  of  a  C-potential  which  is  globally  minimal  (null)  in  the  subset  Cgoai 
of  goal  configurations  and,  except  when  it  tries  to  escape  a  local  minimum,  it  follows  the 
negated  gradient  of  this  C-potential.  On  the  other  hand,  it  uses  a  concept  similar  to  the 
notion  of  valley  for  escaping  encountered  local  minima. 

Starting  from  the  given  initial  configuration,  the  algorithm  follows  the  flow  of  the  negated 
gradient  field  of  the  C-potential  U  until  a  local  minimum,  q/oc,  is  attained  (gradient  motion). 
This  minimum  may  be  a  configuration  where  the  gradient  of  U  is  zero  or  a  configuration 
located  in  the  boundary  of  QCfree  (see  Subsection  4.1).  If  q/oc  £  Cgoah  the  problem  is 
solved.  Otherwise,  the  planner  executes  a  series  of  constrained  motions.  A  constrained 
motion,  denoted  by  M+  (q/oc)  (resp.  Mf  (q/oc)),  with  z  £  [l,n],  consists  of  iteratively  forc¬ 
ing  the  zth  configuration  space  coordinate,  qi ,  to  increase  (resp.  decrease)  by  the  increment 
A i  of  the  grid  QC ,  until  a  saddle  point  of  the  local  minimum  well  is  reached.  Since  q/oc 
is  a  local  minimum  of  the  C-potential,  this  means  that  the  zth  DOF  is  required  to  move 
against  the  force  applied  by  the  C-potential.  At  each  iteration  of  the  constrained  mo¬ 
tion,  the  remaining  n  —  1  coordinates  j  ^  z,  are  selected  according  to  a  best-first  rule. 
Hence,  if  (qq, qi,  ...gn)  is  the  current  configuration,  its  successor  minimizes  the 
C-potential  over  the  set  consisting  of  the  configuration  (gi,  j?n)j  where 

q\  =  qi  +  At-  (if  Mf(qioc)  is  being  executed)  or  qi  -  A*  (if  M“  (q ioc)  is  being  executed)  and 
its  neighbors.  The  motion  thus  tracks  a  kind  of  valley  in  the  (n  —  l)-dimensional  subspace 
orthogonal  to  the  qi  axis.  Hopefully,  it  ultimately  reaches  a  saddle  point  of  the  local  min¬ 
imum  well.  The  planner  recognizes  such  an  event  when  the  C-potential  decreases  again. 
Then,  it  terminates  the  constrained  motion  and  executes  another  gradient  motion. 

Assume  for  an  instant  that  every  constrained  motion  attains  a  saddle  point  of  the  local 
minimum  well.  By  interweaving  gradient  and  constrained  motions,  the  planner  constructs 
a  graph  of  local  minima.  Each  minimum,  after  it  has  been  fully  expanded,  has  2 n  imme¬ 
diate  successors  in  the  graph  and  is  connected  to  each  of  them  by  the  concatenation  of  a 
constrained  motion  and  a  gradient  motion.  Whenever  the  planner  attains  a  minimum,  it 
may  check  if  it  is  a  new  one,  by  comparing  its  coordinates  with  those  of  previously  gen¬ 
erated  minima.  The  search  of  the  local  minima  graph  may  be  conducted  in  a  best-first 
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Figure  20:  Coordinated  Motion  of  4  Mobile  Robots 

fashion  by  iteratively  selecting  a  pending  local  minimum  having  the  smallest  value  of  the 
C-potential  and  generating  its  successors.  In  order  to  avoid  the  systematic  generation  of  all 
the  successors  of  every  selected  pending  minimum,  we  have  implemented  a  variant  of  this 
strategy.  At  each  iteration,  the  implemented  algorithm  selects  a  minimum  q/oc  having  the 
smallest  value  of  U  among  those  whose  successors  have  not  all  been  constructed  yet,  and 
it  generates  a  new  successor  q{oc  of  q ioc.  At  the  next  iteration,  if  U(q{oc)  <  U(q/oc),  the 
algorithm  naturally  selects  q{oc  and  generates  one  of  its  successors;  otherwise,  it  constructs 
another  successor  of  q/oc,  if  any. 

However,  a  constrained  motion  M+  (qjoc)  (or  (q/oc))  may  not  terminate  at  a  saddle  point. 
Instead,  it  may  hit  an  obstacle.  In  that  case,  the  implemented  algorithm  merely  considers 
that  the  local  minimum  has  no  successor  “in  the  direction  of”  increasing  (or  decreasing)  q{. 
An  alternative  would  be  to  treat  the  last  configuration  attained  by  the  constrained  motion 
before  the  collision  as  a  successor  of  q \oc.  This  successor  would  not  be  a  local  minimum, 
but  it  would  nevertheless  be  included  as  a  node  of  the  local  minima  graph  and  then  treated 
as  any  other  node,  except  that  the  constrained  motions  executed  from  it  should  not  be 
commanded  along  the  q{  axis.  It  is  rather  easy  to  construct  examples  where  the  constrained 
motion  technique  succeeds  only  if  these  non-local-minimum  configurations  are  inserted  in 
the  search  graph. 

Without  further  assumptions  on  the  C-potential,  the  above  path  planning  technique  is  not 
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Figure  21:  Coordinated  Motion  of  10  Mobile  Robots 

complete.  Nevertheless,  as  shown  below,  it  produces  interesting  experimental  results. 

7.2  Experimentation 

We  initially  developed  the  constrained  motion  technique  with  the  goal  of  solving  planning 
problems  involving  the  coordinated  motions  of  several  mobile  robots  in  a  workspace  made 
of  narrow  corridors.  Indeed,  consider  the  case  where  a  gradient  motion  leads  two  mobile 
robots  to  roll  in  opposite  directions  in  a  narrow  corridor  where  they  cannot  pass  each  other. 
At  some  point,  they  cannot  progress  further  toward  their  respective  goals  without  colliding. 
Then,  the  gradient  motion  has  reached  a  local  minimum  of  the  C-potential  defined  in  the 
product  configuration  space  of  the  two  robots.  A  constrained  motion  corresponds  to  one  of 
the  robots  moving  backward,  the  other  robot  proceeding  along  the  negated  gradient  of  the 
C-potential,  until  there  is  sufficient  room  for  the  robots  to  pass  each  other. 

We  applied  the  technique  to  plan  the  motion  of  4  robots  in  a  workspace  containing  two 
parallel  walls  forming  a  narrow  corridor  (see  Figure  20).  Each  robot  is  a  L1-disc  of  radius 
a  denoted  by  A{ ,  i  =  which  can  only  translate  in  the  plane.  Its  configuration  is 

defined  as  the  coordinates  (x;,yt)  of  the  center  pi  of  the  square  in  a  workspace  coordinate 
frame15.  The  configuration  space  of  the  whole  robot  system  is  a  8D  space  with  each  con- 

15The  orientation  of  the  coordinate  axes  with  respect  to  the  walls  has  no  major  impact  on  the  operations 
of  the  planner. 
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figuration  represented  by  (xi,  y\, ...,  £4,  j/4).  In  the  example  shown  in  Figure  20,  the  robots 
initially  stand  in  the  corridor  made  by  the  two  walls.  In  the  goal  configuration,  which  is 
uniquely  defined,  the  robots  also  stands  in  the  corridor  but  in  the  reversed  order.  Figure 
20  shows  various  intermediate  configurations  along  the  path  generated  by  the  planner.  In 
this  example,  the  C-potential  was  constructed  as  the  sum  J^t=  1  UPl  (X(pi,  q))  of  the  simple 
W-potentials  computed  at  the  centers  of  the  robots.  The  collision-checking  technique  of 
Subsection  4.1  reduces  to  the  comparison  of  the  L1  distance  di(X(pi, q))  to  the  L 1  radius 
a.  In  the  same  fashion,  the  collisions  between  two  robots  Ai  and  Aj  were  checked  by  com¬ 
paring  the  L1  distance  between  the  points  p,-  and  pj  and  2a.  The  total  planning  time  was 
30  seconds. 

Figure  21  illustrates  another  example  on  which  we  experimented  with  the  planning  tech¬ 
nique.  Ten  mobile  robots,  identical  to  the  robots  of  the  previous  example,  operate  in  a 
workspace  consisting  of  two  rooms  connected  by  corridors.  No  two  robots  can  pass  each 
other  in  any  of  the  corridors,  except  in  the  lateral  corridors  at  the  four  corners,  the  two 
intersections  with  the  central  corridor,  and  the  locations  facing  the  two  doors.  The  whole 
robot  system  has  a  20D  configuration  space.  The  initial  configuration  is  shown  in  the  upper 
lefthand  corner  of  the  figure  and  the  goal  configuration  in  the  lower  righthand  corner.  The 
figure  displays  several  intermediate  configurations  along  the  path  generated  by  the  plan¬ 
ner.  We  first  ran  the  planner  using  the  C-potential  Up,(X(p;,  q)),  as  in  the  previous 
example.  The  planner  was  able  to  construct  a  path,  but  this  path  traversed  several  tens 
of  local  minima.  Many  of  these  minima  were  due  to  the  fact  that  the  robots  were  moving 
in  a  de-organized  fashion  leading  to  intricate  traffic  jams.  Then,  we  ran  the  planner  using 
a  slightly  different  C-potential  obtained  by  adding  to  the  previous  C-potential  a  term  that 
applies  to  each  robot  a  force  perpendicular  to  its  current  direction  of  motion  and  pointing 
towards  its  right.  This  additional  field,  which  is  analogous  to  a  “magnetic  field”,  essentially 
reproduces  the  “drive  on  the  right”  traffic  rule.  In  our  experiments,  it  led  the  various  robots 
to  move  in  a  more  organized  fashion.  Using  this  C-potential,  the  planner  solved  the  problem 
shown  in  Figure  21  by  generating  a  path  traversing  only  a  few  local  minima.  The  total  path 
was  generated  in  about  30  seconds,  which  is  quite  satisfactory  given  the  high  number  of 
robots. 

Since  multi-robot  path  planning  problems  typically  involves  many  DOF’s,  traditional  global 
path  planning  methods  are  often  impractical.  More  specific  approaches  known  as  “prior¬ 
itizing”  and  “velocity  tuning”  have  been  proposed  for  such  problems.  The  prioritizing 
approach  [Erdmann  and  Lozano-Perez,  1986]  consists  of  considering  the  robots  in  sequence. 
When  the  path  of  a  robot  has  been  constructed,  this  robot  becomes  a  mobile  obstacle 
for  the  robots  whose  motions  have  not  yet  been  planned.  The  velocity  tuning  approach 
[Kant  and  Zucker,  1986]  consists  of  first  planning  the  motion  of  each  robot  as  if  it  was 
the  only  robot  in  the  workspace  and  then  tuning  the  velocities  so  that  no  two  robots  col¬ 
lide.  Hence,  both  approaches  first  break  the  planning  problem  for  q  robots  into  q  simpler 
planning  problems,  each  involving  a  single  robot,  and  then  consider  the  interactions  among 
their  solutions.  These  approaches  are  unable  to  solve  multi-robot  problems  involving  tight 
interactions  among  the  robots,  such  as  the  two  problems  shown  above. 

Although  the  constrained  motion  technique  seems  particularly  suitable  for  planning  coordi- 
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nated  motions  of  multiple  mobile  robots,  we  have  also  experimented  with  it  using  the  8-DOF 
manipulator  arm  of  Figure  12  with  the  same  C-potential  as  in  the  experiment  reported  in 
Subsection  5.3.  In  the  example  of  Figure  13,  a  solution  path  was  constructed  in  about  3 
minutes,  which  is  comparable  to  the  result  obtained  with  the  random  motion  technique. 
However,  our  experimentation  with  multi-joint  manipulator  arms  has  been  limited,  and  the 
reliability  of  the  constrained  motion  technique  remains  to  be  verified  for  such  robots. 

8  Conclusion 

In  this  paper,  we  presented  a  collection  of  numerical  potential  field  techniques  for  robot  path 
planning.  All  these  techniques  apply  the  same  general  approach.  They  construct  a  potential 
field  over  the  robot’s  configuration  space,  build  a  graph  connecting  the  local  minima  of  this 
potential,  and  search  this  graph.  The  graph  is  built  incrementally  and  searched  as  it  is 
built.  The  only  precomputation  step  is  aimed  at  the  construction  of  the  local-minimum- 
free  W-potentials  for  the  various  control  points  in  the  robot.  This  step  could  be  avoided, 
but  it  is  an  important  one,  since  it  allows  us  to  construct  “good”  C-potentials.  In  addition, 
this  precomputation  occurs  in  a  space  of  fixed  and  small  dimension  d  =  2  or  3.  In  virtually 
all  practical  cases,  the  size  of  the  graph  of  the  local  minima  is  reasonably  small,  resulting 
in  efficient  path  searching. 

We  proposed  four  techniques  for  constructing  the  local  minima  graph.  These  techniques 
are  based  on  different  ways  of  escaping  the  encountered  local  minima  of  the  C-potential. 
The  best-first  motion  technique  is  very  simple  and  gives  excellent  results  for  robots  with 
few  DOF’s.  However,  it  is  impractical  for  robots  with  many  DOF’s.  The  random  motion 
technique  gives  very  good  results  for  robots  with  few  and  many  DOF’s.  In  addition,  it  is 
highly  parallelizable.  The  valley-guided  motion  technique  gave  significantly  inferior  results. 
However,  it  embeds  ideas  that  one  might  improve  in  the  future  by  exploring  the  concept 
of  valleys  more  thoroughly.  The  constrained  motion  technique,  although  incomplete,  is 
particularly  well-adapted  to  solve  planning  problems  involving  the  coordinated  motions  of 
several  robots.  It  gave  very  good  experimental  results,  specially  for  planning  the  coordinated 
motions  of  several  simple  mobile  robots.  The  principles  underlying  these  four  planning 
techniques  are  somewhat  independent  of  the  representation  of  the  robot’s  workspace  and 
configuration  space.  However,  their  efficient  implementation  was  made  possible  by  the  use 
of  hierarchical  bitmap  representations.  The  planner  implementing  these  techniques  was  able 
to  solve  path  planning  problems,  whose  complexity  measured  by  the  number  of  DOF’s  or 
the  number  of  obstacles  is  far  beyond  the  capabilities  of  previously  implemented  planners. 
For  simpler  problems,  our  planner  is  several  orders  of  magnitude  faster  than  most  previous 
planners. 

We  feel  that  the  techniques  presented  in  this  paper  and  the  experimental  results  obtained 
with  them  make  it  possible  to  realistically  envision  the  development  of  “real-time”  path 
planners.  By  “real-time”,  we  mean  that  the  planners  would  be  able  to  produce  paths  in  a 
very  small  amount  of  time  (say,  a  fraction  of  a  second)  in  almost  all  practical  situations. 
The  construction  of  such  a  real-time  planner  will  probably  require  the  use  of  some  dedicated 
hardware  and  parallel  computing  architecture,  leading  to  the  notion  of  “motion  engine” 
just  like  there  exist  “geometric  engines”  for  performing  graphic  operations  (e.g.,  hidden  line 
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removal).  We  believe  that  the  availability  of  a  real-time  motion  engine  would  open  new 
perspectives  on  some  important  issues  in  industrial  robot  programming  and  autonomous 
robot  control,  and  enable  the  construction  of  efficient  robot  systems  operating  in  partially 
known  and  dynamically  changing  workspaces.  Among  the  four  techniques  presented  in  this 
paper,  the  random  motion  technique  seems  to  provide  the  best  combination  of  qualities 
(time  efficiency,  generality,  reliability).  Since  it  is  also  highly  parallelizable,  it  is  probably 
the  best  candidate  for  developing  such  a  motion  engine.  The  random  motion  technique 
described  in  the  paper  makes  use  of  un-informed  probability  distribution  laws.  If  a  very  fast 
implementation  of  the  technique  was  available,  we  could  envision  the  inclusion  of  a  learning 
module.  In  a  specific  application  domain,  this  module  would  collect  various  pertinent 
statistics,  for  instance  statistics  about  the  distribution  and  the  radii  of  the  minimum  wells, 
and  refine  the  probability  laws  used  by  the  planner,  leading  to  even  faster  and  more  stable 
response  times. 
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